function [Y_datavalues, Base_all, Y_fromeachshocks_all] = zFC_HistoricalDecomposition_partially(Y_used, Autoregressive, B, Resid, StructShocks, t_backuntil, t_until)

k = size(Resid,1);
forhowlong = size(Resid,2);
assert( size(Y_used,2) == size(Resid,2))

[MovAverage_reduced, MovAverage_structural] = zFC_RewriteVAR(Autoregressive, B, forhowlong);

Resid_hd = Resid;
S_hd     = StructShocks;
assert(size(Resid_hd,2) == size(MovAverage_reduced,3))

% t_backuntil = 1;
% t_until     = T_used;

Base_all             = NaN*ones(k,  t_until-t_backuntil+1);
Y_fromeachshocks_all = NaN*ones(k,k,t_until-t_backuntil+1);
Y_datavalues         = NaN*ones(k,t_until-t_backuntil+1);

for step_number = 1:t_until-t_backuntil+1 

    t_at = t_backuntil + step_number - 1; 
    assert(t_at >= t_backuntil) 

    % part of y at time t_at explained by all reduced form shocks occurred between t_backuntil and t_at included
    y_fromallshocks_at = zeros(k,1); 
    for ttt = t_backuntil:t_at
        y_fromallshocks_at = y_fromallshocks_at + MovAverage_reduced(:,:,t_at-ttt + 1)*Resid_hd(:,ttt);

        % checks
        if ttt == t_at
            assert(max(max(MovAverage_reduced(:,:,t_at-ttt+1) - eye(k))) < 0.001)
        end
    end

    % part of y_hat_at explained by each shock considering, for each, all the realizations between t_backuntil and t_at included 
    y_fromeachshocks_at = NaN*ones(k,k);

    for kkk = 1:k
    hd_shock_kkk = zeros(k,1);
      for ttt = t_backuntil:t_at
         hd_shock_kkk = hd_shock_kkk + MovAverage_structural(:,kkk,t_at-ttt + 1)*S_hd(kkk,ttt);

         % checks
         if ttt == t_at
             assert(max(max(MovAverage_structural(:,kkk,t_at-ttt+1) - B(:,kkk))) < 0.001)
         end
      end
    y_fromeachshocks_at(:,kkk) = hd_shock_kkk;

    end
%     assert( max(abs((sum(y_fromeachshocks_at,2) - y_fromallshocks_at))) <     0.0001) % not anymore if some structural shocks are not identified

    % part of y at time t_at not explained by any shock, but by and deterministic part
    base_at = Y_used(:,t_at) - y_fromallshocks_at ;
%     assert( max(abs(( sum(y_fromeachshocks_at,2) + base_at - Y(:,t_at) ))) < 0.0001)

    % store
    Base_all(:,step_number) = base_at;
    Y_fromeachshocks_all(:,:,step_number) = y_fromeachshocks_at;
    Y_datavalues(:,step_number) = Y_used(:,t_at);
    
end
assert( size(Base_all,2) == size(Y_fromeachshocks_all,3) )
assert( size(Base_all,2) == size(Y_datavalues,2) )



end